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ABSTRACT 

The first large-scale quantum mechanical simulations of pulsar glitches are presented, 
using a Gross-Pitaevskii model of the crust-superfluid system in the presence of pin- 
ning. Power-law distributions of simulated glitch sizes are obtained, in accord with 
astronomical observations, with exponents ranging from —0.55 to —1.26. Examples 
of large-scale simulations, containing ~ 200 vortices, reveal that these statistics per- 
sist in the many-vortex limit. Waiting-time distributions are also constructed. These 
and other statistics support the hypothesis that catastrophic unpinning mediated by 
collective vortex motion produces glitches; indeed, such collective events are seen in 
time-lapse movies of superfluid density. Three principal trends are observed. (1) The 
glitch rate scales proportional to the electromagnetic spin-down torque. (2) A strong 
positive correlation is found between the strength of vortex pinning and mean glitch 
size. (3) The spin-down dynamics depend less on the pinning site abundance once the 
latter exceeds one site per vortex, suggesting that unpinned vortices travel a distance 
comparable to the inter-vortex spacing before repinning. 
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1 INTRODUCTION 

Glitches are discrete, randomly timed jumps in the fa- 
mously predictable angular velocity of a pulsar. The statisti- 
cal distributions of glitches in individual pulsars (power-law 
sizes and exponential waiting times) p oint to an underlyin g 
threshold-triggered collective process dMelatos et al.|[ 2008). 
akin t o grain avalanches in sa nd piles ( Wiesenfeld et ahl 



1 19891 : iMorlev fc Schmidt! I1996I) or vorte x avalanches 



type-II superconductors (| Field et all 1 19951 ). In pulsars, the 
'grains' are quantised superfluid vortices. A comprehensive 
theory of why > 10 7 (out of ~ 10 18 ) vortices simultane- 
ously unpin and self-reorganise to create a glitch remains 
elusive. A complete solution relies on connecting the micro- 
scopic (fm-scale) and macroscopic (km-scale) dynamics of 
the neutron superfluid in the pulsar interior. The challenge 
presented by this disparity of scales, especially in replicating 
collective phenomena, has proved formidable. 

In this paper, we present the first comprehensive 
quantum mechanical simulations of pulsar glitches, using 
the Gross-Pitaevskii equation (GPE) to model the tem- 
poral evolution of the superfluid inside a pulsar. The 
GPE is a standard tool in condensed matter physics 
for investigating vortex dynamics in quantum condensates 
ijPethick fc Smithll2002l b Our simulations are unique in four 
ways: (1) the system is much larger than those routinely 
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simu l ated in condensed matter applications (|Tsubota et al.l 
120021: iTsubota et ail l20ld : iPenckwittl I2003L for example), 
and therefore contains many more vortices; (2) in addi- 
tion to the confining potential, we include a grid of pin- 
ning sites to model nuclear lattice defects in the pul- 
sar's crust (|joneslll998bl:lDonati fc Pizzocherol |2004 120061 : 
lAvogadro et all I2007I . |200S| ; iBarranco et alj|20ld ); (3) the 
container and superfluid are coupled via a feedback torque; 
and (4) we allow the angular velocity of the container to 
change non-adiabatically. Thus, the simulations describe 
scales ranging from individual vortex cores to the collec- 
tive behaviour of several hundred vortices in the presence of 
a large grid of pinning sites. The results build on previous 
GPE studie s of the general problem of spas modic superfluid 
spin down (|Warszawski fc Me latos 2010a|) and the micro- 
physics of individual vortex unpinning, including knock-on 
proc esses involving acoustic radiati on and the proximity ef- 
fect (|Warszawski fc Melatoj[2010bl ). 

The paper is structured as follows. Section[2]summarises 
the vortex unpinning paradigm of pulsar glitches, the main 
features of the numerical model, and the feedback mecha- 
nism. In Section [3] we describe our algorithm for extracting 
glitches from simulated spin-down curves. In Section [3] this 
algorithm is employed in a canonical example. In Section [5] 
we construct glitch size and waiting-time distributions from 
relatively small-scale simulations (~ 30 vortices) for differ- 
ent pinning parameters. In Section [6l we vary the inertia of 
the star and the electromagnetic spin-down torque. Section]?] 
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Figure 1. Schematic of how the glitch-finding algorithm oper- 
ates. A and B mark the positions of the two glitches in the inter- 
val 552 < t < 630 on the plot of angular velocity as a function 
of time f!(t). The size of glitch C is the difference in Q between 
B (the local minimum between A and C), and C (the vertical 
distance between the horizontal dashed lines). Simulation param- 
eters: V = 16.6, T) = 1, AVj/Vb = 0.0, R = 12.5, Ax = 0.15, 
At = 5 x 10" 4 , N EM /I C = 10~ 3 , fl(t = 0) = 0.8, n pin /n F = 0.97. 



reports results from larger-scale simulations involving ~ 200 
vortices, in which we look for evidence of collective dynam- 
ics. Parameter studies with the larger-scale simulation are 
too expensive computationally at this stage. In Section[8]we 
synthesise the results and present our conclusions. 



2 VORTEX UNPINNING PARADIGM 
2.1 Overview 

The superfluid in a pulsar is composed of neutron Cooper 
pairs, which i nterpenetrate the nuclear lattice in the crust 
(Migdal ll959| y The superfluid attempts to match the rota- 
tion of the crust by forming quantised vortices, each carrying 
a quantum of circulation n = h/(2m), whose areal density 
is proportional to the angular velocity of the superfluid. In 
the absence of pinning, vortices form a hexagonal Abrikosov 
lattice, causing the superfluid to rotate as a rigid body. In 
equilibrium, the mean areal density of vortices is np = 20/ K. 
In reality, however, vortices pin at or between nuclear lattice 
sites and/or defects, thus distorting the equilibrium vortex 
configuration. The extent of distortion depends on the con- 
figuration and abundance of pinning centres supplied by the 
crustal lattice of nuclei. Pinning prevents the superfluid from 
decelerating smoothly with the crust. Glitches are believed 
to occur when a large number of vortices simultaneously 
unpin and move outward t o new pinned positions. Indeed , 
GPE simulations reported in|Warszawski & Mel atos! (|2010al ) 
exhibit stick-slip vortex dynamics and resulting glitch-like 
interruptions to the smooth spin down of the rotating con- 
tainer. Similar spin-up events were observed in laboratory 
experiments using rotating c ontainers filled with helium II 
jTsakadze fe Tsakadzd[l980h . 



We construct a simple, two-dimensional pulsar model 
comprising a superfluid inside a rotating, decelerating crust, 
in the presence of a pinning grid that co-rotates with the 
crust. We include self-consistent coupling between the su- 
perfluid and crust by ensuring that angular momentum is 
conserved overall. 

A popular explanation for pulsar glitches is that an- 
gular momentum is transferred from the superfluid to the 
crust in response to sudd en unpinning, outward motion, 
and repinning of vortices (|Packardl Il972l ; lAnderson fc Itohl 
1975). Crust-superfluid differential rotation produces a fic- 
titious Magnus force, Fm, that unpins vortices when it ex- 
ceeds some threshold. Fm is transverse to the relative veloc- 
ity between the vortex line and superfluid flow; it is similar 
in origin to the Bernoulli lift force in conventional hydrody- 
namics. For a vortex line moving with velocity v through a 
superfluid of density p and velocity v s , the Magnus force is 
given by 



F M = p« x (v - v s 



(1) 



where n is directed along the vortex axis. For a pinned 
vortex, v is the crust velocity; for a vortex at r = b, one 
has v = Q.bi. For N v vortices in an Abrikosov lattice, the 
Feynman relation approximates (the result is exact in the 
infinite- vortex limit) the azimuthal component of the super- 
fluid velocity at radius r as v$ ~ nN v /(2irr). 

The vortex unpinning paradigm presumes that a cat- 
alyst exists to spontaneously unpin between 10 7 and 10 15 
vortices during a single glitch, with the range deduced from 
measurements of glitch sizes. To date, the catalyst has not 
been identified securely. If vortices unpin independently, the 
number of unpi nning events in any t ime interval follows Pois- 
son statistics (|Hakonen et al.lll99ST ). On its own, a Poisson 
process does not account for the observed power-law distri- 
bution of glitch sizes. Hence, a mechanism is needed which 
correlates vortices, so that the unpinning rate of a particular 
vortex increases when another vortex (either nearby or far 
away) unpins; collective unpinning avalanches are essential 
to explaining pulsar g li tches. In recent GPE simulations, 
Warszawski & Melatos (201Q3) demonstrated two mecha- 
nisms - in addition to the growing global velocity shear as 
the crust spins down - for triggering and sustaining vortex 
avalanches: (i) acoustic radiation from moving vortices, and 
(ii) vortex proximity knock-on, as unpinned vortices move 
past their pinned neighbours. 

Acoustic radiation res ults from th e propagation of 
Kelvin waves on vortices JVinenl l200ll ), from vortex re- 
connection jLeadbeater et al]|200lh and due to dissipation 
by ac celerating vortices ( Parker et al.|[20o3 : iBarenghi et al.l 
l2005h . The disturbance in the density field as the sound 
waves propagate is capable of unpinning weakly pinned vor- 
tices. The amplitude of the disturbance is a function of vor- 
tex acceleration. Therefore, the strongest acoustic radiation 
results from the repi nning phase, as a vortex spira ls in to 
the new pinning site jWarszawski fc Melatosll2010bh . When 
the level of dissipation in the system is high, the acoustic 
disturbance generated by a repinning vortex is weaker than 
for a low-dissipation system, and hence it is harder to unpin 
vortices acoustically. 

The second effect, proximity knock-on, occurs when 
two vortices approach sufficiently closely that the Mag- 
nus force on at least one of them exceeds the unpinning 
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Quantity 


Simulation 


Simulation 


Pulsar 


Pulsar 




(dimensionless) 


(dimensional) 


(dimensionless) 


(dimensional) 


R 


12.75 


10- 16 m 


10 21 


10 4 m 


V 


16.6 


10 4 MeV 


10" 3 


10° MeV 


n 


0.8 


10 24 Hz 


10" 24 - 10~ 21 


10" 1 - 10 2 Hz 


Nem/Ic 


10" 3 


10 45 Hz s- 1 


10-63 _ 1Q-59 


10" 15 - 10- 10 Hz s- 1 


*total 


10 3 


10~ 21 s 


10 32 - 10 33 


10 s - 10 9 s 


lip 


0.13 


2.8 x 10- 16 m- 2 


10 s - 10 12 


10~ 8 - 10~ 4 m" 2 


^pin 


0.19 


2.6 X 10- 16 m- 2 


10° - 10 4 


10-14 - 10-1° m- 2 


N v 


30 - 200 


30 - 200 


10 16 - 10 19 


10 16 - 10 19 


V 


1 


1 


10~ 6 


10~ 6 


Ax 


0.15 


1.5 x 10~ 16 m 






At 


5 x 10" 4 


5 x 10" 28 m 







Table 1. Dimensionless and dimensional parameters for a typical simulation (first and second columns) and pulsar (third and fourth 
columns). Dimensionless quantities are scaled so as to r ecover Eq. ^ from the dimensional form of the GPE, making the replacements 
described in Section 2 of Warszawski & Mclatos (2010a). Pulsar quantities are quoted as a range representative of the pulsar population. 
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Figure 2. Timing signature of a simulated glitch. Top: Angular velocity of the crust as a function of time Q(t) for f2(t = 0) = 0.8. Q(t) is 
smoothed with a top- hat function of width At sm = 1.0. The arrows mark the positions of glitches identified using the automated glitch 
finder described in Section [3] Bottom left: Close-up of a glitch at t = 736, bracketed by the dotted rectangle in the top panel. The dotted 
curve represents the linear model fitted to the pre-glitch spin-down curve. Bottom right: Phase residual after removing a pre-glitch linear 
fit to the interval 710 < t < 730 [U(t = 710) = 0.942, tl = -8.98 X 10~ 4 ] from the interval 710 < t < 770 shown in the bottom left panel. 
Simulation parameters: Vq = 16.6, r? = 1, AVi/Vo = 0.0, R = 12.5, Ax = 0.15, At = 5 X 10~ 4 , N EM /Ic = 10" 3 , n pin /n F = 0.97. 
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Figure 3. Superfluid angular momentum as a function of time (L z )(t), in units of h/2, for the canonical spin-down experiment described 
in Section [4] The arrows indicate the position of glitches identified using the glitch-finding algorithm defined in Section [3] Simulation 
parameters: t] = 1, AVi/V = 0.0, R = 12.5, Ax = 0.15, At = 5 x 10" 4 , N EM /I C = 10" 3 , fi(t = 0) = 0.8, n pin /n F = 0.97. 



threshold. It was found that a ballistic unpinned vortex 
can unpin a pinned vortex at a larger radius via the prox- 
imity effect. The distance of closest approach between the 
two vortices is inversely prop ortional to pinning strength 
l|Warszawski fc Melatoil2010bh . 



2.2 Gross-Pitaevskii model 

In this section we summaris e the quantum mechanical pulsa r 
model first introduced in IWarszawski fc Melatosl (|2010al ). 
The interested reader who wishes to reproduce our results 
is referred to the latter paper for technical details. In the 
model, a two-dimensional, zero-temperature Bose-Einstein 
condensate represents an equatorial projection of the stellar 
superfluid. The time evolution of the superfluid density \ip\ 2 , 
defined in terms of the order parameter ip(x.,t), in a poten- 
tial V, with chemical potential a, observed in a reference 
frame rotating with angular velocity Q, is described by the 
dimensionless GPE 



(i 



l) d JL = 
y) dt 



QL z tp 



(2) 



In arriving at Eq. (2j, we define the sound speed c 3 — 
(nop/m) 1 / 2 , the healing length (also the characteristic 



length-scale) £ = h/ {mnog) 1 ^ 2 , the characteristic time-scale 
h/(nog), and the energy scale nog, where g is the self- 
interaction strength, m is the mass of each boson (twice the 
neutron mass), and no is the mean particle density. The an- 
gular momentum operator is given by i z = —izd/d(j> (with 
units h/2), where <f> is the azimuthal angle. V is the sum of 
the confining potential representing the crust and the pin- 
ning potential. 



In line with standard practice (I Jin et al.l 



Choietal. 199; 



Kasamatsu ct al 



Tsubota et~aIll2002l ; [Penckwitt et all 



1996; 



2002; 



we include in Eq. ([2]) a phenomeno- 



logical treatment of dissipation, modelled by the term 
—jdtfj/dt. Its effect is to suppress sound waves, thus re- 
ducing the relaxation time-scale of the superfluid and en- 
suring that the wave function can adjust to non-adiabatic 
ch anges in Q, due to glitch e s. On t he other hand, as reported 
in I Warszawski fc Melatosl (|2010bl ). sound waves play an im- 
portant role in vortex unpinning avalanches, by facilitating 
knock on. Therefore, one must be careful not to raise 7 to the 
point where it quells the latter effect. An accurate estimate 
of 7 in a pulsar is not available at present. 

The dissipationless GPE (7 = 0) models a superfluid at 
T = 0, which does not contain viscous excited-state com- 
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Figure 4. Angular velocity as a function of time Q(t) for the canonical example described in Section^ smoothed with a top-hat window 
function of width At sm = 0.0, 1.0, 5.0 and 10 (solid, dotted, dashed, dot-dashed curves respectively; the curves are shifted vertically 
down by 0.0, 0.1, 0.2 and 0.3 respectively to improve visibility). Arrows below each curve mark the positions of glitches found using 
the glitch-finding algorithm described in Section [3] as At sm increases, the number of found glitches decreases. Simulation parameters: 
Vo = 16.6, ri = 1, AVi/Vo = 0.0, R = 12.5, Ax = 0.15, At = 5 x 10~ 4 , JV EM // C = 10~ 3 , n(t = 0) = 0.8, n pin /n F = 0.97. 



ponents. In contrast, solutions to Eq. J5J f or 7^0 are by 
definition at non-zero temperature. IPenckwitt et alj (|2002l ) 
attributed dissipation in experiments to atom transfer be- 
tween the viscous thermal cloud and the condensate. Exper - 
iments with Bose-Einstein condensates (|Halian et al.ll200lh 
showed that the thermal cloud also influences vortex nucle- 
ation. A self-consistent description o f this phenomenon re- 
quire s additional terms in the GPE l|Wouters fc Carusottol 
12007ft . which are not included in most studies, including 
ours. Laboratory experi ments claim 7 = 0.03 in helium II 
l|Abo-Shaeer et ail[20 02 r ); we use 7 = 0.05. The dissipative 
term drives particle loss from the system, which we correct 
by adjusting the chemical potenti al, such that the num- 
ber of particles remains cons tant (|Kasamatsu et all 120031 ; 
IWarszawski fc Melatodl2010al ). 



2.3 Container and pinning potential 

Using a fourth-order Runge-Kutta algorithm in time, and a 
fourth-order finite difference scheme for the spatial deriva- 
tives, we solve Eq. ([2]) for a superfluid in a rotating potential 
V , in the frame co-rotating with the potential. 



V is the sum of a confining potential, representing the 
edge of the pulsar inner crust, and the pinning potential, 
Vpin, representing the nuclear lattice in the crust, which co- 
rotates with the confining potential. V p i n comprises a square 
grid of equidistant spikes, each one described by 

V pta ,i(r) =Vi[l + tanh(A|r-Ri|)] , (3) 

where R; is the position of the centre of the pinning site, 
and Vi and A -1 parametrise its strength and width respec- 
tively. In all cases, the spikes in the pinning grid are equally 
spaced, with a mean areal density n p i n . In Section [5.31 we 
use the ratio n p i n /np as a measure of the relative abun- 
dance of pinning sites and vortices. In Section f5. 5 1 we draw 
Vi from a top-hat distribution of pinning strengths with 
mean Vo and width AVi\ in all other simulations, Vi is 
constant (Vi = Vo). A non-zero range of pinning strengths 
(AVi/Vo ~ 1) is likely to play an important role in glitch 
dynamics. Each small unpinning event biases the distribu- 
tion of p inning strengths of occupied pinning sites towards 
large Vi (|Newman fc Sneppenll 19961 ; iMelatos fc Warszawskil 
120091 ). until a large global velocity shear or strong acoustic 
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Figure 5. Glitch statistics as a function of the degree of smoothing for the spin-down experiment in Fig. [4] Left: Probability density 
function (pdf) of fractional glitch sizes p(Afi/fi), binned logarithmically and graphed on log- log axes, for top- hat window functions of 
width At sm = 0.0, 1.0, 5.0 and 10 (solid, dotted, dashed, dot-dashed curves respectively). Right: Cumulative distribution of glitch waiting 
times P(At) for top-hat window functions of width At sm = 0.0, 1.0, 5.0 and 10 (solid, dotted, dashed, dot-dashed curves respectively), 
graphed on linear axes. Simulation parameters: Vq = 16.6, rj = 1, AVi/Vb = 0.0, R = 12.5, Ax = 0.15, At = 5 X 10 — 4 , N^m/Ic = 10 -3 , 
fl(t = 0) = 0.8, ra pin /n F = 0.97. 



perturbation comes along and unpins many strongly pinned 
vortices simultaneously. 

The parameters of vortex pinning with in a pulsar have 
been studied theoretically in great detail ([Blasio &: Lazzaril 
199cj;|jonej|l998bl;[Hir"asawa fc Shibazakj|200j ; |jonedl2003[ 
Avoeadro et all 120081 : iBarranco et all |201oF l Of particu- 



to 



lar interest is the pinning strength, the density of pin- 
ning sites, and the region within the star where pin- 
ning is strongest (the superfluid pairing state changes with 
depth). The strength of pinning is usually calculated from 
the energy difference between a vortex sitting atop a col- 
umn of spherical defects and a vortex positioned inter- 
stitially. In general, vortices pin to nuclear sites to min- 
imise the region from which superfluid is excluded (the vor- 
tex core is devoid of superfluid). Depending on the super- 
fluid density, however, pinning can also occur inters t itially 



nuid density, however, p inning can als o occur interstitially 
(iPizzochero et all \l99T\: iDonati fc Pizzocherol |2004 I2006J ; 
lAvogadro et al.1 [20071 ; iPizzocherol 120071 ). For a pulsar, the 
characteristic size of the pinning site (the width of a nu- 
cleus) is approximately equal to the superfluid healing length 
(~ 1 fm), which is significantly smaller than the lattice 
(~ 10 fm) and inter-vortex (~ 10 12 fm) spacings. 



2.4 Feedback torque 

Vortex motion changes the superfluid angular momentum 
(L z ), the expectation value of the operator L z . The response 
of the crust to these changes is an essential element of any 
pulsar glitch theory. We employ a simple conservation ar- 
gument, wherein changes in (L z ) are communicated instan- 
taneously to the crust. In reality, changes in the superfluid 
flow are communicated t o the crust b y sound and Kelvin 
waves along vortex lines (|Joneslll998ah , but these processes 
are fast. Given an electromagnetic (EM) spin-down torque, 
TVem, the angular velocity of the crust, Q, evolves according 



Jc 



d(L z ) 



N E 



(4) 



dt dt 
where 7 C is the moment of inertia of the crust, and is 
the same angular velocity that appears in the GPE. We 
emphasise that, even for smooth deceleration, Nem/Ic does 
not equal the observed pulsar spin-down rate, as this rate 
is the combination of two effects: the external EM torque, 
and the spin-up re sponse to gradual (non-glitch) spin down 
of the superfluid (Warszawski & Mcla toa l2010al ). We also 
note that if the vortex lattice is significantly distorted by 
pinning, the superfluid does not rotate rigidly. Hence it is 
not meaningful to assign to it a unique angular velocity or 
moment of inertia. We solve Eq. @ simultaneously with 
Eq. @ in the following sections. 



3 AUTOMATED GLITCH FINDER 

Our first task is to define exactly what we mean by a glitch. 
In practice, glitches in pulsar data are extracted 'by hand', 
by looking for discontinuities i n the slope and/or cu r vature 
of the pulse phase. Recently, IChukwude &: Uramal (|2010T i 
employed an automated glitch-finding algorithm, based on 
a least-squares goodness-of-fit parameter, to identify micro- 
glitches. Since our simulations provide extremely high time 
resolution (QAt ~ 10~ 4 revolutions/time step compared to 
£lAt ~ 10 5 revolutions/time step for the observations), we 
resolve the glitch spin-up phase; we can identify glitches di- 
rectly by looking for spin-up events, rather than resorting 
to polynomial fitting of the phase residuals. The difference 
between simulation and astrophysical time scales is also re- 
flected in which physical processes can be resolved. In our 
simulations, we track the motion of individual vortices be- 
tween pinning sites, whereas observed glitches appear as un- 
resolved spin-up events, which likely result from a multi- 
generation vortex cascade. 
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Figure 7. Vortex configurations and spin-down curves for different pinning strengths Vq. Top: Greyscale plots of the final (t = 1000) 
superfluid density for V; = 0, 8.3, 16.6 and 33.7 (left to right respectively). The density greyscale runs from dark (low density) to light 
(high density). Bottom: Angular velocity of the crust as a function of time Q(t) for Vi = 0, 8.3, 16.6 and 33.7 (solid, dotted, dashed 
and dot-dashed respectively). Simulation parameters: rj = 1, AVi/Vo = 0.0, R = 12.5, Ax = 0.15, At = 5 X 10 -4 , (Vem// c = 10 — 3 , 
fl(t = 0) = 0.8, n pin /n F = 0.97, At sul = 1.0. 



For the purposes of this investigation, we employ the 
following simple algorithm for identifying glitches in fi(t) 
curves. 

(1) Smooth £l(t) with a top-hat window function w(t) of 
width Atsm, viz. f2 sm (t) = f* dt'w(t - t')fi(t'). 

(2) Calculate Af2 S m(ii) = fism(ii) — f2 sm (t;_i) for all dis- 
crete time points ti in the simulated range. 

(3) Find all t s for which Af2 sm (i s ) > 0. For each t s , 
find the smallest t g (> t a ) such that AO sm (i 9 ) ^ and 
Af2 sm (t 9 +i) < (points A and C in Fig. Q] meet this cri- 
terion), and record Q aul (t g ). 

(4) Locate the global f2 sm (t) minimum between i 9 _i and 



t g , occurring at i g , m m (point B in Fig.[T}. Define glitch size as 
AQ/£l = [Q sl n(t g ) — Q S m(t gtin i n )]/Q(tg) (the vertical distance 
between the dashed lines in Fig. [T}, and waiting time as 



At = tg 

and C in Fig. HJ . 



Ig-X 



(the horizontal distance between points A 



The principal weakness of this algorithm is the stipulation in 
step (3) that Af2 sm > 0, which does not capture decreases 
in (L z ) that are smaller than A^EMAt. Depending on the 
moment of inertia of the crust, AS7 g due to real glitches 
may not be greater than the negative Af2 sm due to the spin- 
down torque [(Nem/ Ic)&t], and hence this algorithm is not 
suitable for observational data. 
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Figure 8. Glitch statistics for different pinning strengths corresponding to the f2 curves graphed in Fig. [7] Left: Probability density 
function (pdf) of fractional glitch sizes p(Af2/f2) for Vb = 0, 8.3, 16.6 and 33.7 (solid, dotted, dashed and dot-dashed respectively) 
logarithmically binned and graphed on log-log axes. The dotted curve is a power-law least squares fit to the dotted histogram with 
index —0.818. Right: Cumulative probability of glitch waiting times P(At) graphed on log-linear axes for Vb = 0, 8.3, 16.6 and 33.7 
(solid, dotted, dashed and dot-dashed respectively). Simulation parameters: rj = 1, AVi/Vo = 0.0, R = 12.5, Ax = 0.15, At = 5 X 10 -4 , 
Nem/Ic = 10~ 3 , Q(t = 0) = 0.8, n pin /n F = 0.97, At sm = 1.0. 
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Figure 6. Close-up of glitches identified using the glitch-finding 
algorithm described in Section [3] as a function of the degree of 
smoothing. The plots cover the interval 356 < t < 386 (taken from 
the Q curve graphed in Fig. |3}. The curves are smoothed with a 
top-hat window function of width At Bm = 0.0, 1.0, 5.0 and 10 
(solid, dotted, dashed, dot-dashed curves respectively). The num- 
ber and position of glitches, shown as arrows in the top panel, 
changes with At sm ; for Ar. sm ^ 5, one glitch is detected. Simu- 
lation parameters: V = 16.6, 77 = 1, AVi/V = 0.0, R = 12.5, 
Ax = 0.15, At = 5 x 10" 4 , N EM /I C = 10" 3 , Q(t = 0) = 0.8, 
n P in/"F = 0.97. 



The arrows under the spin-down curve in the top panel 
of Fig. [2] indicate the positions of glitches found using the 
above algorithm for a curve smoothed with a top-hat win- 
dow function of width At sul — 0.1. The glitch sizes and 
waiting times are then used to construct probability density 
functions. In Section [4] we discuss the two bottom panels of 
Fig. [2] which demonstrate how we subtract a linear fit to 



the pre-glitch spin-down curve (bottom left panel), and the 
resulting phase residual curve (bottom right panel). 



4 CANONICAL GLITCH 

Let us now examine a canonical example, which demon- 
strates the ability of the model in Section [2] to generate as- 
trophysically relevant results. Table [1] lists the set of canon- 
ical parameters for this investigation. In the first column, 
the simulation parameters are given in dimensionless units; 
their dimensional counterparts are provided in the second 
column. The third and fourth columns quote the range of 
values associated with a typical pulsar in dimensionless and 
dimensional units respectively. 



4.1 Set up 

The nuclear lattice of the neutron star inner crust is repre- 
sented by an 11 x 11 square grid of equal-strength pinning 
sites (Vb = 16.6, AVi/Vb = 0.0), within a box of side length 
15, with crust radius R = 12.5 (in dimensionless units). The 
simulation grid is a 200 x 200 rectangular grid; the greyscale 
plots appear circular because l^ 2 = outside the crust. The 
time step and grid spacing are At = 5 x 10 -4 and Ax — 0.15 
respectively. 

The initial conditions are established by accelerating 
the crust from rest to Q = 0.8 instantaneously and allowing 
the system to reach a steady state (defined as when the 
fractional change in total energy during one time step drops 
below 10 -5 ). Once a steady-state solution is obtained, the 
external spin-down torque N-em{~ 10 _3 / c ) is switched on, 
and the sup erfluid- crust system is allowed to evolve. 
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62.31 


0.334 
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Table 2. Summary of simulation results for different values of the smoothing width At sm [corresponding Q(t) curves shown in Fig. [4]. 
N g is the number of glitches detected, (AQ/Q) is the mean glitch size, (At) is the mean waiting time, and r p is the Pearson correlation 
coefficient; r p = 1, —1 and indicate positive, negative and no correlation respectively. A is the mean glitch rate that parametrises the 
cumulative distribution of waiting times, 7 is the power-law index that parametrises the pdf of glitch sizes. The K-S probability pxs i s 
given in parenthesis; the fit can be rejected at the 1 — Pks confidence level. Simulation parameters: Vq = 16.6, 7] = 1, AVi/Vo = 0.0, 
R = 12.5, Ax = 0.15, At = 5 x 10~ 4 , N EM /I C = 10~ 3 , fl(t = 0) = 0.8, n pin /n F = 0.97. 



4.2 Spin-down curve 

The spin down observed after Nem switches on is spas- 
modic. By watching movies of how the density field \ip( x , t)\ 2 
evolves, we can easily see why: pinning prevents vortices 
from moving smoothly towards the boundary as the crust 
decelerates; instead, vortices hop between pinning sites, re- 
sulting in step- like decreases in (L z ). A graph of (L z ) versus 
time is presented in Fig. [3] The positions of glitches identi- 
fied using the automated glitch-finding algorithm described 
in Section [3] are marked with arrows. Equation tells us 
that each step down in (L z ) is accompanied by a step up in 
O (so long as the change in (L z ) exceeds N-^yiAi). We see 
this clearly by comparing the top panel of Fig. [2] and Fig. [3] 

A close-up of one of the glitches (bracketed by the dotted 
rectangle in the top panel of Fig. [2j is shown in the bottom 
left panel of Fig. [2] Does this look like a pulsar glitch? An 
important difference is that the spin-down rate before and 
after the glitch does not change; the Q.(t) curve after the 
glitch parallels the linear timing model describing the pre- 
glitch spin down. This contrasts with some (not all) observed 
glitches, which freq uently produce a p ermanent increase in 
the spin-down rate (|Wong et al.ll200ll ). The simulation out- 
put also lacks a post-glitch recovery phase, during which 
spin down is faster than the global average. Given that time 
scales associated with the post-glitch recover y phase, i.e. 
~ 10 29 h/(nog), [typically lasting days to weeks (|Wong et al.l 
2001)] are associated with coupling of the superfluid to a 
viscous fluid component n ot included in these simula tions, 
this result is unsurprising l|van Evsden fe Melatos![2010l ). In 
our simulations, the time for Q(t) to return to its pre- 
glitch value is proportional to the size of the glitch and 
is comparable to the sound-c rossing of the system, R/c B 
jWarszawski fc Melatoll2010al ). 

Practical pulsar glitch detection involves comparing 
measured pulse times of arrival (TOAs) against TOAs pre- 
dicted by a spin-down model. When a glitch occurs, pulses 
arrive more frequently than the model predicts, and the 
phase residuals (the difference between the number of rev- 
olutions predicted by the timing model and the number of 
revolutions inferred from the TOAs) turn negative and grow 
linearly with time. In the bottom right panel of Fig. [2l we 
plot the difference between the simulation output and a tim- 
ing model defined by 

fimodclW = 0(* = 0) +tlt , (5) 

for the interval bracketed by the dotted rectangle in the top 
panel. The timing model, with slope Q = —8.98 x 10 -4 , has 



been extracted by performing a least squares fit of Eq. 
to Q(t) (710 ^ t ^ 730) and is graphed as a dotted curve in 
the bottom left panel of Fig. [2] The phase residual Acj>, in 
units of radians, is therefore 

A0(t)= f Sl modtA (t)dt - f Q(t)dt . (6) 
Jo Jo 

The shape of the residual curve, graphed in the bottom 
right panel of Fig. [2l qualitatively resembles observations; 
see Fig. 1 in lShemar fc Lvnel (1996), for example. 

4.3 Statistics 

We now construct probability density functions (pdfs) of 
sizes and waiting times using the automated glitch finder. 
Motivated by the hypothesis that glitches are produced 
by avalanche dynamics, we look for evidence that the 
pdfs of sizes and waiting times are distributed as power- 
laws [p(Afi g /fi) oc (Afig/fi) 7 ] and exponentials [p(Ai) = 
Aexp(— AAt)] respectively. Evidence is buildi ng that pdfs of 
this f orm yield good fits to observational data jMelatos et al.l 
2008). To fit these functional forms, we find the parameter 
that minimises the D statistic (the maximum distance be- 
tween the cumulative distribution of the model and that 
of the data). The fit can be rejected with 1 — pks confi- 
dence, where pks is the Kolmogorov-Smirnov probability 
jPress et al.ll 19921 ). 

To begin with, we note that the smoothing width At sln 
changes the glitch statistics, as summarised in Table [2] 
Glitch positions are marked on the spin-down curve in Fig. [4] 
using arrows for At sul = 0.0, 1.0, 5.0 and 10.0 (solid, dotted, 
dashed and dot-dashed curves respectively). The curves have 
been separated vertically for clarity. As At sm increases, the 
number of glitches decreases from N s = 105 for At snl — 0.0 
to Ng = 13 for At snl = 10.0. The mean waiting time (At) 
and glitch size (Afi/fl) increase five- and nine- fold respec- 
tively. The trend in (Af2/fi) versus At sln is evident in the 
size pdfs graphed on log-log axes in the left panel of Fig. [5] 
the pdfs for longer At sln are weighted further towards the 
large-(Afi/n) end. Notably, the best-fit power-law index, 7 
(we warn the reader that this is not the same quantity as 
the dissipation parameter 7 that appears in Eq. [2|, does 
not change dramatically with At sm ; we find 7 = —0.994 for 
Atsm = (7 can be rejected at the 25.5% confidence level), 
and 7 = —1.104 for At sul = 10.0 (rejected at the 5.3% con- 
fidence level). 

On the other hand, increasing At sln changes markedly 
the shape of the cumulative waiting-time distribution (right 
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panel of Fig. [5] graphed on linear axes). For At sm ^ 1.0, 
the cumulative probability is dominated by short waiting 
times (i.e. convex down). For Atsm ^ 5.0, waiting times are 
more evenly distributed. Exponential fits to the cumulative 
probability return mean glitch rates between A = 0.011 (for 
At sm = 10.0; rejected at the 87% confidence level) and A = 
0.87 (for At sm = 0.0; rejected at the with near certainty). An 
exponential functional form for the cumulative distribution 
for At sm = 1.0 is rejected with the least confidence (55.1%). 

Glitch positions and distributions depend on At sm be- 
cause smoothing fi(t) affects the morphology of individual 
glitches. Understanding this tendency is essential to inter- 
preting observational data properly using our algorithm. For 
example, increasing At sm reduces the size of a glitch, and 
in some cases erases the glitch altogether. To visualise this, 
in Fig. [6] we graph an interval of Q(t) containing an obvious 
bump (when viewed on the scale graphed in Fig- [4j for differ- 
ent smoothing widths. The arrows in the top panel mark the 
positions of glitches identified by the automated algorithm 
for different At sm . In the unsmoothed (At sm = 0.0, solid) 
and Atsm — 1-0 (dotted) curves, we see two obvious bumps 
at t — 371 and 367.5 (for At sm = 0.0 the algorithm finds 
8 glitches that are not visible by eye on the scale plotted). 
However, for At sm = 5.0 and 10.0, the two bumps merge, 
and the algorithm finds only one glitch. For At sm = 10.0, 
the glitch is ~ 16% smaller than for At sm = 5.0. Increased 
waiting times are a direct corollary of glitches being erased 
by smoothing. 

After watching movies of \4>\ 2 (x,t), we conclude that 
glitches too small to be discerned by eye on the scale de- 
picted in Fig. [6] are associated with 'jiggling', not unpin- 
ning. By 'jiggling' we mean vortices wobble around their 
pinning site without unpinning. Multi-peaked glitches arise 
from multi-stage vortex motion (i.e. one vortex moves, caus- 
ing a second vortex to move a short time later), precisely 
the collective avalanche behaviour we are seeking. Observa- 
tionally, we can only detect pulsar glitches resulting from 
the motion of many vortices (the smallest observed glitches 
correspond to ~ 10 7 unpinnings). For this reason, we hence- 
forth adopt Atsm = 1.0, which also happens to produce 
the most exponential P(At). Thus we retain multi-peaked 
glitches like the one in Fig. [51 whilst excluding glitches not 
associated with unpinning. 

An important recent result from pulsar data analysis is 
the absence of a r eservoir effect: gl i tch size is uncorrela ted 
with waiting time l|Wong et al.ll200ll ; [Meiatos et alj200Sl ). In 
Table [2] we list the Pearson correlation coefficient, r p , relat- 
ing An/fi and At as a function of At sm . For At sm < 5.0, the 
two quantities are uncorrelated (r p = —0.057 and 0.050 for 
At sm = 0.0 and 1.0 respectively). For At sm 5.0, r p differs 
significantly from zero (0.428 and 0.334 for At sm = 5.0 and 
10.0 respectively), suggesting that sizes and waiting times 
are correlated. The interpretation of this correlation depends 
on the glitch definition. Since the correlation emerges for 
smoothing scales that erase multi-peaked glitches and 'jig- 
gling', it is likely that this correlation exists if the glitch- 
finding algorithm does not distinguish between single- and 
multiple-vortex motion. Simulations large enough to house 
several isolated regions of unpinning activity should not ex- 
hibit such a correlation dMelatos et al.l l2008). 

As observational duty cycles and measurement preci- 
sion improve, e.g. with the advent of multibeam telescopes, 



astronomers will begin to resolve the spin-up phase of a 
glitch, multi-peaked glitch es, and post-glitch oscillations 
l|Kramer fc Stappersi l2O10h . Our results suggest that im- 
proved measurement precision (equivalent to reducing At sm ) 
should lead to a large increase in the number of observed 
glitches, without altering the shape of the size distribution. 
For example, reducing At am five-fold doubles the number 
of glitches. Conversely, the mean glitch rate increases. More 
speculatively, we also predict that the low-Af2/fi end of the 
glitch size pdf is dominated by ongoing, non- collective vor- 
tex creep fsee lAlpar et al.lll984l . for example), likely charac- 
terised by a Gaussian cut off. 



5 PINNING PHYSICS 



Systematic differences in the pinning strength from pul- 
sar to pulsa r may result from different cooling histories 
(I Jones! 120011 ). which alter the morphology and abundance 
of crystal defects. The role of the pinning strength Vo in 
glitch physics is crucial. If pinning is weak, vortices spread 
out smoothly and homologously, allowing the superfluid to 
decelerate in sympathy with the crust; differential rotation 
does not build up in stress reservoirs. On the other hand, 
if pinning is strong, it may stiffen the vortex lattice to the 
point where the crust cracks before the vortices are unpinned 
by the Magnus force (|Horowitz fc Kadau|[2009h . Intuitively, 
we expect that stronger pinning results in larger glitches; ei- 
ther vortices move further once they have unpinned, or more 
vortices unpin in order to reverse the crust-superfluid shear. 
In this section we use simulations to test these intuitive re- 
lations between glitch distributions and pinning strength. 

5.1 Equal pinning potentials 

We begin by examining how spin down depends on Vo, in 
the situation where all the pinning sites are of equal strength 
(AVi — 0). The four greyscale plots of superfluid density at 
the top of Fig. [7] are the final states (t — 1025) for sim- 
ulations with Vo — 0, 8.3, 16.6 and 33.3 (left to right re- 
spectively). Vortices are easily identified as darker spots. 
The lighter spots represent unoccupied pinning sites. All 
vortices, except for one in the north-east quadrant of the 
second panel, are pinned. 

Without pinning, the vortex lattice expands smoothly 
and homologously, (L z ) decreases smoothly, and the super- 
fluid exerts a continuous spin-up torque on the crust via 
Eq. @. The spin- up torque lessens the gradient of (the 
black curve in Fig. [7]) when compared with no- feedback, lin- 
ear spin down (the solid grey curve). Small bumps in the 
black curve in Fig. [7] are numerical noise. 

We now switch on pinning, via an 11 x 11 equidistant, 
equal-strength grid. (L z ) exhibits step- like decreases, accom- 
panied by jumps in Q,, in the dotted, dashed and dot-dashed 
curves in Fig. (7) The first glitch (at ti) occurs later when 
the pinning is stronger, viz. ti = 200 for Vo = 8.3 compared 
to ti = 440 for Vo = 33.3. The delay can be understood by 
considering the velocity shear between pinned vortices and 
the superfluid in the absence of feedback: 

Av(r, t) = «s - [«(t = 0) - (N EM /Ic)t]r . (7) 
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Table 3. Glitch size and waiting-time statistics for different pinning potentials Vo corresponding Q(t) curves shown in Fig. [7] N s is 
the number of glitches detected, (Afi/O) is the mean glitch size, (At) is the mean waiting time, and A is the activity parameter. A 
is the mean glitch rate that parametrises the cumulative distribution of waiting times, 7 is the power-law index that parametrises the 
pdf of glitch sizes. The K-S probability pks is given in parenthesis; the fit can be rejected at the 1 — pks confidence level. Simulation 
parameters: n = 1, AVJ/Vfc = 0.0, R = 12.5, Ax = 0.15, At = 5 x 10" 4 , N EM /I C = 10" 3 , Sl(t = 0) = 0.8, n pin /n F = 0.97, At sm = 1.0. 



Whilst vortices remain pinned, v s is approximately constant. 
Therefore, it takes longer for the Magnus force to unpin 
strongly pinned vortices. 

Even in the strongest pinning case, before any unpin- 
ning has occurred, the spin-down curve is less steep than 
Q(t) = tN-Ehi/Ic (the solid grey curv e in Fig. P. The reason 
is subt le; it is discussed in detail in IWarszawski fc Melatoi 
l|2010al ). While remaining pinned, vortices respond to in- 
creasing differential rotation by migrating smoothly to the 
outer edge of their pinning site. Hence the crust experiences 
a gradual spin-up torque. Put another way, pinning is not 
a binary condition; it is weakened gradually as the shear 
grows. 

At the end of the simulation, S7 is highest (fl = 0.29) for 
the weakest pinning and lowest (fi = 0.16) for the strongest 
pinning, because stronger pinning potentials can sustain 
greater le vels of stress. [Note that this is not the reser- 
voir e ffect l|Link et al.ll 19991 ; I Wong et alJbOOlljMelatos et all 
2008), which predicts a correlation between glitch size and 
waiting time not found in observational data.] In other 
words, the areal density of remaining vortices after a fixed 
amount of time is proportional to the pinning strength. 

At this point we draw attention to vast discrepancies 
between typical pulsar parameters and the physical scales 
modelled by our simulations. The extremely short simula- 
tion time ( 10 -21 s in physical units), and the extremely 
high angular velocity (fl = 10 24 Hz in pulsar units com- 
pared to the typical pulsar angular velocity S7 ~ 10 Hz) 
represent highly non-adiabatic conditions. 

5.2 Statistics 

A key goal of this paper is to understand the observed scale 
invariance of glitch sizes, w hich is believed to a rise from 
collective vortex behaviour (Me latos et al.l I2008T ). Global 
properties of the glitch distribution, such as the number of 
glitches, mean glitch size, and mean waiting time for each 
Vo studied in Section 15.11 are listed in Table [3] along with 
the activity parameter 

A = 7V g ^ASl g /^At , (8) 

used by pulsar astronomers to characterise th e total glitch- 
induced spin up of a pulsar (|Lvne et al .11 19981 ). 

In Fig. [8] we graph logarithmically binned size pdfs on 
log-log axes. Our glitch-finding algorithm identifies 16, 48, 
23 and 26 glitches for Vo = 0, 8.3, 16.6 and 33.3 respectively. 
The pdfs for simulations with pinning (Vo > 0) exhibit a 
bias towards larger glitches for larger Vo, which also shows 



up as a ~four-fold increase in mean glitch size for Vo = 
33.3 compared to Vo = 8.3. Power-law fits to the pdfs yield 
similar indices (listed with the K-S probability in parenthesis 
in the penultimate column of Table [3| for 8.3 ^ Vo ^ 33.3. 
However, the confidence level with which the fits can be 
rejected ranges from 23.8% for V = 16.6 to 70.7% for V = 
33.3. We remind the reader that glitches detected in the 
Vo = case arise from numerical noise, rather than vortex 
unpinning; this assertion is supported by the small mean 
glitch size (10 4 (AS7/S1) = 13.00). 

In the right panel of Fig. [8] we graph cumulative distri- 
butions of waiting times for the simulations shown in Fig. 
Mean glitch rates, A, derived from e xponential fits , are lis ted 
in the final column of Table [3] In iMelatos et all (|2008l ). it 
was found that waiting times for seven of the nine pulsars 
that have glitched more than six times are well represented 
by exponential distributions. Taking the results in Fig. [8] all 
together, we suggest that A does not depend monotonically 
on Vo because there is another degree of freedom, (Afl/fi), 
through which the additional stress can be released. 

The activity parameter does not appear to scale mono- 
tonically with Vo. However, A almost doubles between Vo = 
8.3 and Vo = 33.3. We do not have a physical explanation 
for this behaviour at this time. 

5.3 Pinning site abundance 

If pinning occurs at individual nuclei, then ~ 10 9 pin- 
ning sites separates neighbouring vortices in a pulsar. Al- 
ternatively, if pinning occurs at grain boundaries or along 
macroscopic faults in the crustal lattice, the ratio of pinning 
sites to vortices is not necessarily large, and pinning sites 
may not be evenly spaced. In this section we explore the 
role of the pinning site density (n p i n sites per unit area) rel- 
ative to vortex density (riF = 2S7/k vortices per unit area). 
We present results for simulations with: (i) few or no pin- 
ning sites, (ii) n p i n ~ tlf, and (iii) many more pinning sites 
than vortices. We parametrise the pinning site abundance 
in terms of the ratio n p i n /np in the initial state. 

Four pinning scenarios are illustrated in the top panels 
of Fig. [9] as greyscale plots of \ip\ 2 (x,t = 1000) for n p i n /nF 
increasing from left to right. The middle panels both cor- 
respond to regime (ii). In every panel except the leftmost 
one (n p in/nF = 0.11), all vortices in the final state are 
pinned. For n p i n /riF = 0.11, even in the final state, there 
are only enough pinning sites to pin nine out of 17 vor- 
tices (i.e. there are no unoccupied pinning sites). Notably, 
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Table 4. Glitch size and waiting-time statistics for different n p ; n /np [corresponding Q(t) curves shown in Fig. 0. N g is the number of 
glitches detected, (AO/O) is the mean glitch size, (At) is the mean waiting time, and A is the activity parameter. A is the mean glitch 
rate that parametrises the cumulative distribution of waiting times, 7 is the power-law index that parametrises the pdf of glitch sizes. 
The K-S probability pxs is given in parenthesis; the fit can be rejected at the 1 — Pks confidence level. Simulation parameters: Vj = 16.6, 
T] = 1, AVi/Vo = 0.0, R = 12.5, Ax = 0.15, At = 5 x 10" 4 , N EM /I C = 10" 3 , Q(t = 0) = 0.8, Ai sm = 1.0. 



for npi n /nF = 0.11, the vortex lattice is regularly spaced, 
suggesting that pinning alters the orientation, but not the 
geometry, of the vortex lattice. For n p i n /np = 0.52 and 
Wpin/nF = 0.97, vortices are irregularly spaced, as they are 
forced to adhere to the geometry of the pinning grid. When 
pinning sites significantly outnumber vortices (right panel, 
npm/?iF = 1.90), the vortex lattice is once again evenly 
spaced, because enough pinning sites are available for a 
vortex to 'choose' a site near its equilibrium (Abrikosov) 
position. It should be noted that this equilibrium position 
changes as the crust spins down. 

The bottom panel of Fig. [9] graphs the spin-down curves 
for the four pinning abundances illustrated in the top panels. 
The example with few pinning sites (n P i n /nF = 0.11, solid 
curve) is relatively smooth. Although our algorithm detects 
25 glitches, they are small ((Afi/n) = 14.84 x 10~ 4 ; cf. 
(Afi/fi) = 13.00 x 10~ 4 for Vo = 0) and well separated 
((At) — 37.10), suggesting that they arise in the rare event 
that an outwardly moving vortex runs over a pinning site. 
By contrast, for n p i n ~ tif (dotted, dashed and dash- dotted 
curves), the spin down is roughly linear until vortices begin 
to unpin at t ~ 275, whereupon vortices move towards the 
edge in spasmodic jumps between pinning sites (vortices do 
not stop at all pinning sites in their path). For all n p i n /riF > 
0.11, the gradual outward migration of vortices across their 
pinning sites reduces Q, even before vortices have unpinned. 
For n p i n /riF = 0.11, f2 is considerably less than Nem/Ic for 
t > 50, as the unpinned vortices move smoothly outward 
almost as soon as the spin-down torque is switched on. 

It is physically interesting that the total spin down 
of the crust over long time periods is almost identical for 
w p in/nF > 1 and n p i n /nF ~ 1. The aggregate effect of 
glitches is insensitive to increases in n p i n /nF over and above 

unity. One interpretation is to say that the average dis- 

1/2 

tance travelled by a vortex once it unpins is of order n F , 
rather than n pi ]/ 2 , so that vortices maintain an approximate 
Abrikosov lattice, without developing the large-scale spatial 
inhomog eneities required in self-organi sed critical avalanche 
models (|Warszawski fc M elatos 2008). However, this hy- 
pothesis is challenged by three trends. First, the number of 
glitches and activity parameter increase for n p i n /n,F > 0.11; 
these qualities scale with n~^J 2 , even though rip 1 ^ 2 is ap- 
proximately fixed. Second, the mean waiting time is in- 
versely proportional to n p i n /nF; see Table [4] for details and 
the right panel of Fig. [9] for cumulative distributions of At. 
Third, the mean glitch size does not change monotonically 
with n p i n /nF. Taken together, these results suggest that the 
typical distance travelled by a vortex in a glitch is a play off 



between the number of vortices that unpin, the distance to 
the next available pinning site, and the distance the vortex 
must move to reverse the shear that unpinned it. 

5.4 Maximum size 

Glitch size is governed by the total change in (L z ) and hence 
depends on the number of vortices that unpin, the distance 
that each vortex travels before repinning, and the location 
from which it unpinned. From the pdfs of Afi/fl in Fig. 1101 
we see that the largest glitch is approximately the same for 
all ^pin ^ np. If vortex motion is determined by the first 
encounter with a pinning site, rather than the imperative to 
reduce shear and maintain equal spacing, then we expect the 
mean glitch size to scale inversely with the density of pinning 
sites. This is not observed. By contrast, maximum glitch size 
does scale monotonically with Vb [max(10 4 Afi/fl) = 36, 538, 
1057, and 1641 for Vb = 0, 8.3, 16.6 and 33.3 respectively], 
because vortices need to move further to reduce the shear 
sustained by stronger pinning. 

5.5 Unequal pinning potentials 

Pinned vortices at a given radius experience an equal Mag- 
nus force. If every pinning site is of equal strength, then 
vortices at a given radius should unpin simultaneously, yield- 
ing periodic glitches of equal size, contrary to obs ervations 
|Melatos et alJl200Sl : fMelatos fc Warszawskj|2009l ). In this 
section, we investigate the effect of varying the width of the 
pinning strength distribution AVi/Vo. The glitch statistics 
from these simulations are summarised in Table [S] 

The spin-down curves in Fig. [11] demonstrate that vary- 
ing AVi/Vo makes little difference to the total spin down 
over the time scale of the simulation; Q(t — 1000) is ap- 
proximately the same for all ^ AVi/Vo ^ 1. However, the 
timing of the first departure from linear spin-down scales 
as t g i « 1250 + 200(1 - AVi/Vo). In addition, the glitch 
statistics change systematically with AVi/Vo. The pdfs of 
Afl/fi graphed in the left panel of Fig. 1121 confirm that the 
minimum glitch size depends on AVi/Vo, via the minimum 
pinning strength. The correlation is not strict because the 
number of draws from the pinning strength distribution is 
relatively small (11 x 11 = 121), so the minimum V, does 
not monotonically scale with AVi/Vo. We detect no system- 
atic change in the exponent of power-law fits to the size pdf. 
However, we note that the fits can be rejected with high con- 
fidence for AVi/Vo = (78.8%) and AVi/Vo = 1.0 (74.7%). 
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Figure 9. Results from numerical experiments for different pinning abundances n p ; n /riF. Top: Greyscale plots of the final superfiuid 
density for n p i n /rip = 0.11, 0.52, 0.97 and 1.90 (left to right respectively). The density scale runs from dark (low density) to light 
(high density). The feint dots are unoccupied pinning sites; darker dots are pinned vortices. Bottom: Angular velocity of the crust as a 
function of time f2(t) for n p ; n /np = 0.11, 0.52, 0.97 and 1.90 (solid, dotted, dashed and dot-dashed respectively). Simulation parameters: 
Vi = 16.6, rj = 1, AVi/Vo = 0.0, R = 12.5, Ax = 0.15, At = 5 x 10~ 4 , N EM /Ic = 10~ 3 , fl(t = 0) = 0.8, At sln = 1.0. 



Cumulative waiting time distributions are graphed in 
the right panel Fig. [12] Exponential fits reveal a tendency 
for the mean glitch rate to increase with AVi/Vo (A = 0.07 
for AVi/Vo = 1.0 and 0.043 for AVi/Vo = 0.0). This find- 
ing disagrees with the uniform pinning scenarios discussed 
in Section 15.11 and Section 15.21 in which A does not scale 
inversely with Vo- 

To interpret the above findings, we suggest that, even 
for AVi/Vo = 0, the effective potential at each pinning site 
is a function of its departure from an equilibrium position in 
the Abrikosov lattice. If a vortex is slightly inside or outside 
(radially) its Feynman position, the effective pinning energy 



is decreased or increased respectively, and hence the distri- 
bution of pinning strengths can be treated as heterogeneous, 
even when AV5 = 0. In a pulsar (n p i n /fiF ~ 10 14 if all nu- 
clear sites can pin; see TablefJ}, the pinned vortices resemble 
an Abrikosov array. Hence we expect that AVi/Vo is more 
influential in a pulsar than in the simulations presented here. 

6 STELLAR PARAMETERS 
6.1 Crust moment of inertia 
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Figure 10. Glitch distributions for spin-down experiments with different pinning abundances npi n /np for the H curves graphed in 
Fig. [9] Left: Probability density function (pdf) of fractional glitch sizes p(Afi/fi) for n p i n /np = 0.11, 0.52, 0.97 and 1.90 {solid, dotted, 
dashed and dot-dashed respectively) logarithmically binned and graphed on log-log axes. The solid grey curve is a power-law least squares 
fit to the solid histogram with index —0.801. Right: Cumulative probability of glitch waiting times p(Ai) graphed on log-linear axes 
for ripin/np = 0.11, 0.52, 0.97 and 1.90 (solid, dotted, dashed and dot-dashed respectively). Simulation parameters: Vq = 16.6, 7] = 1, 
AVi/Vo = 0.0, R = 12.5, Ax = 0.15, At = 5 x 10~ 4 , N EM /I C = 10~ s , Q(t = 0) = 0.8, At sm = 1.0. 
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N s 


10 4 (Afi/T2) 


(At) 


A 


7(pks) 


A(PKS) 


0.0 


28 


231.2 


21.72 


59.82 


-0.818(0.212) 


0.048(0.360) 


0.25 


28 


173.3 


22.13 


44.13 


-0.676(0.596) 


0.043(0.712) 


0.5 


27 


231.0 


26.74 


47.19 


-0.805(0.504) 


0.038(0.232) 


0.75 


21 


261.8 


33.42 


37.80 


-0.776(0.671) 


0.051(0.008) 


1.0 


30 


206.2 


26.45 


47.10 


-0.901(0.253) 


0.070(0.055) 



Table 5. Glitch size and waiting-time statistics for different ranges of pinning strengths AVi/Vo [corresponding f!(t) curves shown in 
Fig. 1111 . jVg is the number of glitches detected, (AQ/Q) is the mean glitch size, (At) is the mean waiting time, and A is the activity 
parameter. A is the mean glitch rate that parametrises the cumulative distribution of waiting times, 7 is the power-law index that 
parametrises the pdf of glitch sizes. The K-S probability pxs ' s given in parenthesis; the fit can be rejected at the 1 — pxs confidence 
level. Simulation parameters: V = 16.6, 77 = 1, n pin /n F , R = 12.5, Ax = 0.15, At = 5 X 10~ 4 , Af EM // c = 10~ 3 , Q(t = 0) = 0.8, 
At sm = 1.0. 



The inertia of the pulsar crust, relative to the inertia 
of the pinned superfiuid, controls the crust's response to 
changes in (L z ). Let 77 — I c /Is be the ratio of the crust's 
moment of inertia, I c , to the moment of inertia that the su- 
perfiuid would have if it were a rigid body. If the crust is 
light (small 77), small decreases in (L z ) produce large crustal 
spin-up events. For sufficiently sma ll 77, even vortex creep 
jAlpar et all Il984l : iLink et allll993t ) should produce small 
but observable glitches. For larger 77, on the other hand, 
crustal spin up is smalle r . In w hat follows, we adopt the view 
proffered bv lLink et ail |l999Tl that the 'crust' is composed 
of the rigid nuclear lattice, all charged fluid components that 
are coupled electromagnetically to the rigid lattice, and the 
unpinned superfiuid (which may be a large fraction of the to- 
tal). Hydrodynamic models o f glitch recovery suggest 77 ~ 1 
(van Evsden fc Melatosll2010F ). Results for the numerical ex- 
periments presented in this section are summarised in Ta- 
bled 

As 77 increases, the mean glitch size diminishes; a given 
change in (L z ) results in a smaller change in Q. For example, 
from the spin-down curve in Fig. 1131 we see that doubling 
77 halves Q(t = 1000). Put another way, as 77 increases, un- 



pinning and outward vortex motion do progressively less to 
reduce the crust-superfiuid shear. Hence, the timing of the 
unpinning events changes with 77. Cumulative waiting time 
distributions are graphed in the right panel of Fig. 1141 They 
show that the mean waiting time decreases as 77 increases. 
The waiting-time distribution is weighted to the low At end 
for 77 = 2.0 compared to 77 = 1.0 (e.g. (At) = 27.93, 25.19 
and A = 0.035, 0.065 for 77 = 1.0, 2.0). 

The resistance of a heavy crust to spin up also produces 
smaller glitches. This result is confirmed by considering the 
pdf of Afi/fl for two cases, 77 = 1.0 and 2.0. Looking at the 
left panel in Fig. 1141 we find that the heavier crust is accom- 
panied by smaller glitches, e.g 10 4 (Af2/fi) = 188, 74.73 for 
77=1.0,2.0. 

The vertical distance between the curves in Fig. [13] in- 
creases with time, demonstrating that the increased glitch 
rate is not sufficient to compensate for reduced glitch sizes. 
Over longer time periods than those simulated here, Cl self- 
regulates so that the crust-superfiuid shear exactly equals 
the value required to sustain the average unpinning rate 
leading to Cl (given Nem and 77). The activity parameter 
A provides a means of comparing the relative changes in 
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Figure 11. Angular velocity of the crust as a function of time f2(t) for spin-down experiments with AVi/Vo = 0.0, 0.25, 0.5, 0.75 and 1.0 
(solid, dotted, dashed, dot-dashed and triple- dot- dashed respectively). Simulation parameters: Vo = 16.6, r\ = 1, AVi/Vo = 0.0, R = 12.5, 
Ax = 0.15, At = 5 x 10" 4 , N EM /I C = 10" 3 , fi(t = 0) = 0.8, Ai sm = 1.0. 




Figure 12. Glitch distributions for spin-down experiments with different AVi/Vo for the Q(t) curves graphed in Fig. 1111 Left: Probability 
density function (pdf) of fractional glitch sizes (AQ/Q) for AVi/Vo = 0.0, 0.25, 0.5, 0.75 and 1.0 (solid, dotted, dashed, dot-dashed and 
triple- dot- dashed respectively) logarithmically binned and graphed on log-log axes. The solid grey curve is a power-law least squares 
fit to the solid histogram with index —0.818. Right: Cumulative probability of glitch waiting times p(At) graphed on log-linear axes 
for AVi/Vo = 0.0, 0.25, 0.5, 0.75 and 1.0 (solid, dotted, dashed, dot-dashed and triple- dot- dashed respectively). Simulation parameters: 
Vo = 16.6, n = l, R = 12.5, Ax = 0.15, At = 5 x 10" 4 , N EM /I C = 10" 3 , Q(t = 0) = 0.8, At am = 1.0. 
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Figure 13. Angular velocity of the crust as a function of time Q(t) for spin-down experiments with lighter (7) = 1) and heavier (r] = 2) 
crusts (solid and dotted curves respectively). Simulation parameters: Vq = 16.6, AVi/Vo = 0.0, -R = 12.5, Ax = 0.15, At = 5 X 10 — 4 , 
Nem/Ic = 10~ 3 , 0(1 = 0) = 0.8, At sm = 1.0. 




Figure 14. Glitch distributions for spin-down experiments with different AVi/Vo for the Q(t) curves graphed in Fig. 1131 Left: Probability 
density function (pdf) of fractional glitch sizes (AO/fi) for r\ = 1 and 2 (solid and dotted curves respectively) logarithmically binned 
and graphed on log-log axes. The dotted grey curve is a power-law least squares fit to the dotted histogram with index —0.846. Right: 
Cumulative pdf of glitch waiting times p(At) graphed on log-linear axes for rj = 1 and 2 (solid and dotted curves respectively). Simulation 
parameters: V, = 16.6, AVi/Vo = 0.0, R = 12.5, Ax = 0.15, At = 5 x 10" 4 , JVem / ic = 10" 3 , Q(t = 0) = 0.8, At sm = 1.0. 
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JV g 


10 4 (AQ/n) 


(At) 


A 


7(PKS) 


A(pks) 


1.0 


27 


188.0 


27.93 


36.45 


-0.846(0.615) 


0.035 (0.450) 


2.0 


37 


74.73 


25.19 


22.00 


-0.760(0.916) 


0.065 (0.152) 



Table 6. Glitch size and waiting-time statistics for different crustal moments of inertia r\ [corresponding Q(t) curves shown in Fig. 113] . 
N g is the number of glitches detected, (Af!/f2) is the mean glitch size, (At) is the mean waiting time, and A is the activity parameter. 
A is the mean glitch rate that parametrises the cumulative distribution of waiting times, 7 is the power-law index that parametrises the 
pdf of glitch sizes. The K-S probability pxs ' s given in parenthesis; the fit can be rejected at the 1 — pxs confidence level. Simulation 
parameters: V = 16.6, AVi/V = 0.0, n pin /n F , R = 12.5, Ax = 0.15, At = 5 x 10~ 4 , N EM /I C = 10~ 3 , U(t = 0) = 0.8, Ai sm = 1.0. 



response to the spin-down torque in (At) and (AQ/Q); A 
decreases by ~ 33% as 77 doubles from r\ — 1.0 to 2.0, demon- 
strating again that the increased glitch rate does not com- 
pensate for the large decrease in glitch size. 

6.2 Electromagnetic spin-down torque 



The magnetic dipole torque acting on a pulsar ulti- 
mately drives the crust-superfiuid shear essential to the 
(un)pinning model of pulsar glitches. Here, we quantify 
the torque in terms of an electromagnetic spin-down rate, 
Nem/Ic, in the absence of feedback. The larger the torque, 
the more quickly Magnus stresses build up in the pinned 
vortex lattice. Results from simulations in which Nem/Ic is 
varied are summarised in Table [7] 

In Fig. [15] we graph spin-down curves for Nem /Ic = 
10~ 2 ' 5 (solid), 10~ 3 (dashed), 10~ 3 ' 5 (dotted) and 10~ 4 
(dot-dash). For Nem/Ic = 10 _4 '°, we detect few glitches 
(Ng = 3), owing to the small shear that accumulates over 
the total lifetime of the simulation [£total-^EM/-Tc = 0.4]. 
Importantly, the first glitch does not appear at the same 
value of Q for different spin-down torques; Q(ti) = 0.45, 
0.46, 0.61 and 0.44 for Nem/Ic = lO" 4 , 10~ 3 ' 5 , 10~ 3 
and 10 -2 ' 5 respectively. This result reinforces t he fin dings 
of hysteresis reported in Jackson fe Barenghil l|2006l ) and 
IWarszawski fc Melatosl (|2010al ): the spin-down history, not 
merely the current Q, determines the number and position 
of vortices. 

For a crust-superfiuid system with r\ — 1.0, the ob- 
served crust spin-down rate is Aem/(27c); see Eq. Q 
with dLz/dt = I c dQ./dt. For N s 3> 1, the results in 
Fig. [TS] demonstrate that the fraction of the torque- 
driven spin down t tota \NEM/Ic that is reversed by glitches 
decreases as Nem/Ic increases. Quantitatively, we find 
I c {AQ g )/((At)NEM) = 0.269, 0.245 and 0.137 for Nem/Ic = 
10~ 3 ' 5 , 10~ 30 and 10" 4 respectively. We interpret this to 
mean that weaker torques drive the system adiabatically, al- 
lowing the superfluid to match the crust's deceleration. In a 



pulsar, the slow spin down (0,/Q < — 10~ s - ) corresponds 
to the Iow-TVem regime. Hence, we can expect on large time 
scales the superfluid deceleration matches that of the crust, 
thus avoiding unfettered growth of differential motion. 

Size and waiting-time distributions for Nem/Ic > 



10' 



are plotted as logarithmically binned pdfs and cumu- 



lative probability distributions in the left and right panels 
of Fig. 1161 respectively. As Nem/Ic increases, one obtains 
more frequent, smaller glitches. The results reported in Ta- 
ble [7] show that the slopes of power-law fits to the size pdf 
also increase monotonically with Nem/Ic, as does the rate 



extracted from exponential fits to waiting-time pdfs (from 
A = 0.0289 to 0.112 for Nem/Ic = 10~ 3 ' 5 and Nem/Ic = 
10 -2 ' 5 respectively; since the Nem/Ic = 10 -4 curve contains 
only 3 glitches, we do not trust the fitted A). Systematic in- 
creases in the activity parameter A (from A = 0.0958 to 
73.98 for Nem/Ic = 10~ 40 and Nem/Ic = 10~ 2 ' 5 respec- 
tively) confirm that increasing Nem/Ic does not trigger the 
same vortex movements as for smaller Nem/Ic- 



7 LARGE-SCALE SIMULATIONS 



The scale-invariant statistics of pulsar glitches are ex- 
pected to emerge most clearly in the many- vortex limit. In 
this section we increase the system size, and hence the vor- 
tex number, by an order of magnitude. Prevented by com- 
putational expense from undertaking a full parameter study, 
we present examples of two systems with uniform and vari- 
able pinning strengths (AVi/Vo — 0.0 and 0.8 respectively), 
each with ~ 200 vortices in total. The simulations ran for 
~ 2 weeks on a single 2.5 GHz processor. Statistics from the 
simulations are summarised in Table [8] 

To maintain the same spatial resolution as in the sim- 
ulation, we increase the crust radius to R = 38 and impose 
a 66 x 66 pinning grid (with Vo = 36). We also increase 
the number of particles in the system 10- fold (via no), to 
provide sufficient contrast to distinguish vortices from the 
background density. We initialise the experiment by instan- 
taneously accelerating the crust from rest to SI = 0.8 and 
letting the system reach a steady state. Once a steady state 
is achieved, the crust is subjected to a spin-down torque 
Nem = 10~ 3 J C . 

The greyscale plot in the left panel of Fig. [T7] depicts 
the initial superfluid density for uniform pinning (AVi/Vo = 
0.0). The two insets zoom in on regions where vortices sit 
interstitially. The smooth increase in \ip\ 2 with radius, in- 
dicated by the radial colour gradient, represents the stan- 
dard parabolic density profi le, analogous to the meniscus o f 
a rotating bucket of water (Warszawski & Mela tosl bOlOal ). 
Initially the vortices (dark dots in Fig. 1171 light dots are 
unoccupied pinning sites) are surrounded by a vortex-free 
corona in the region 4 ^ r ^ 38, highlighting the distor- 
tion and hy steresis wrought by the pinnin g grid (see Sec- 
tion IVB of I Warszawski fc Melatol l2010al ). By the end of 
the simulation, the vortices are approximately evenly dis- 
tributed across the pinning grid; see the greyscale plot in 
the right panel of Fig. 1171 The radial density gradient is 
shallower because Q, is lower than its initial value. We can 
see that vortices have moved from their initial positions; 
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Table 7. Glitch size and waiting-time statistics for different Nem/Ic for the Q(t) curves shown in Fig. 1151 N g is the number of glitches 
detected, (Af2/f2) is the mean glitch size, (Ai) is the mean waiting time, and A is the activity parameter. A is the mean glitch rate that 
parametrises the cumulative distribution of waiting times, 7 is the power-law index that parametrises the pdf of glitch sizes. The K-S 
probability pxs i s given in parenthesis; the fit can be rejected at the 1 — Pks confidence level. Simulation parameters: Vq = 16.6, 77 = 1, 
AVi/Vo = 0.0, n pin /n F , R = 12.5, Ax = 0.15, At = 5 x HT 4 , AT EM // C = 10" 3 , Q(t = 0) = 0.8, Ai sm = 1.0. 




Figure 15. Angular velocity of the crust as a function of time Q(t) for spin-down experiments with different electromagnetic spin-down 
torques Nem/Ic = — 10 -2 ' 5 , — 10 — 3 '°, — 10 -3 ' 5 and — 10~ 40 (solid, dotted, dashed and dot-dashed curves respectively). Simulation 
parameters: V = 16.6, q = 1, AVj/Vb = 0.0, R = 12.5, Ax = 0.15, Ai = 5 x 10~ 4 , fl(t = 0) = 0.8, At sm = 1.0. 



vortices present in the insets in the left panel of Fig. \T7\ are 
absent from the right panel. 

The spin-down curves for AVi/Vo = 0.0 and 0.8 (solid 
and dotted curves in Fig. 1181 respectively), smoothed with 
a top-hat window function of width At sm = 0.1, reveal 
only small glitches (10 4 (Af2 g /ng) = 0.752 and 0.591 for 
AVi/Vo — 0.0 and 0.8 respectively). This result is unsurpris- 
ing. In larger systems, each vortex corresponds to a smaller 
fraction of the total angular momentum. To achieve glitches 
as large as those produced by a R = 12.5 system, ~ 10 



times more vortices must unpin simultaneously. Dissipation, 
which suppresses sound waves, renders knock-on events less 
likely; hence such large-scale collective unpinning is rare in 
the systems simulated here. The insets show close-ups of the 
largest glitch (top right and bottom left for AVi/Vo = 0.0 and 
AVi/Vo = 0.8 respectively). The spin- up phase of the largest 
glitch in the AVi/Vo = 0.0 curve is bumpy. In fact, the ver- 
tical distance between the bumps is of order the height of 
the largest glitch in the AVi/Vo = 0.8 curve, which does not 
exhibit the same 'bumpiness'. 
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Figure 16. Glitch distributions for spin-down experiments with different electromagnetic spin-down torques Wem/^c for the f2(t) curves 
graphed in Fig.[l5] Left: Probability density function (pdf) of fractional glitch sizes (Afi/O) for Nem/Ic = — 10 -2 ' 5 , — 10 -30 , — 10 -3 ' 5 
and — 10~ 40 (solid, dotted, dashed and dot-dashed curves respectively) logarithmically binned and graphed on log-log axes. The solid 
grey curve is a power-law least squares fit to the solid histogram with index —0.984. Right: Cumulative probability of glitch waiting 
times p(At) graphed on log-linear axes for Nem/Ic = — 10 — 2 ' 5 , — 10~ 30 , — 10 -3 ' 5 and — 10 — 4 (solid, dotted, dashed and dot-dashed 
curves respectively). Simulation parameters: V = 16.6, g = 1, AVi/V = 0.0, R = 12.5, Ax = 0.15, At = 5 X 10~ 4 , U(t = 0) = 0.8, 
Af sm = 1.0. 
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Table 8. Glitch size and waiting-time statistics for different At sm and AVj/VJj [corresponding Q(t) curves shown in Fig. 1181 for large-scale 
simulations containing ~ 200 vortices. 7V g is the number of glitches detected, (AQ/Q) is the mean glitch size, (At) is the mean waiting 
time, and A is the activity parameter. A is the mean glitch rate that parametrises the cumulative distribution of waiting times, 7 is 
the power-law index that parametrises the pdf of glitch sizes. The K-S probability pks is given in parenthesis; the fit can be rejected 
at the 1 — Pks confidence level. Simulation parameters: Vq = 36, g = 1, AVi/Vo = 0.0, n p i n /np, R = 38, Ax = 0.15, At = 5 X 10 — 4 , 

n E m/Ic = 10- 3 , n(t = 0) = 0.8. 



Pdfs of Af2/S! and cumulative distributions of At are 
graphed in the left and right columns of Fig. Ir respectively; 
the solid curves are for AVi/Vo = 0.0 and the dotted curves 
are for AVi/Vo = 0.8. Since there is an obvious turnover in 
the size distribution at Afl/Q ~ 7 x 10~ , power-law fits 
are to p(Af2/f2 ^ 7 x 10~ 7 ). The top panels refer to glitches 
detected in the unsmoothed Q curve, whereas the bottom 
panels are smoothed with At sm = 0.1. For both values of 
AVi/Vo, smoothing flattens the size distribution. Compari- 
son with Fig. [5] reveals that the smallest glitch in the larger 
system is six times smaller than the smallest glitches for 
the R = 12.5 systems discussed in previous sections, which 
agrees with the ~six— fold increase in vortex number N v . If 
a single vortex is responsible for the smallest glitch, then 
the size of the smallest glitch should scale inversely with 
N v . For At sm = 0.1, the fitted power law is marginally shal- 
lower for AVi/Vo = 0.8 (7 = -1.05) than for AVi/Vo = 0.0 
(7 = -1.23). 

The waiting times tend to lengthen as AVi/Vo in- 
creases, as observed in smaller systems (see Fig. 1121) . For 
AVi/Vo = 0.0, smoothing broadens the waiting-time distri- 
bution significantly. The top right panel of Fig. [12] reveals 



a strong periodicity at At ~ 0.05 in the waiting-time dis- 
tribution for At sm = 0.0. This may be evidence of the pe- 
riodic glitches expected from a system in which stress in- 
creases over time by the same amount on vortices pinned 
with equal strength. Despite the overall slower spin down 
for AVi/Vo = 0.8 (cf. Fig. 1180. the longer waiting times are 
not matched by larger glitches. This is likely due to one of 
two reasons: (i) much of the crust-superfluid shear is reversed 
by smooth vortex motion, unhindered by weak pinning; or 
(ii) glitches are too small to cause a bump (positive change) 
in Q, and hence are not detected by our glitch-finding algo- 
rithm. 

Two significant differences between this experiment and 
those described in Section 2] — Section [5] prevent us from at- 
tributing the altered spin-down behaviour (smaller glitches) 
entirely to system size. First, the ratio of pinning sites to 
vortices is significantly higher (n P m/np » 5) in this ex- 
periment than in the smaller systems, whilst the spacing 
of the pinning sites, relative to the vortex core size, is 
smaller. It may be that such a closely spaced pinning grid 
presents the vortices with a nearly homogeneous, elevated 
potential, effectively weakening the pinning strength. Sec- 
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Q = 0.800 



Q = 0.200 





Figure 17. Greyscale plots of the initial (Q = 0.8, left) and final (f2 = 0.2, right) superfluid density |i/>| 2 for a 600 X 600 simulation with 
R = 38 and Vi = 36. The insets zoom in on regions with dimensions 3x3. Vortices are recognisable as the darker dots, initially clustered 
towards the centre of the crust and evenly spaced at the end of the simulation. Unpinned vortices can be seen in the initial but not the 
final insets. Simulation parameters: V = 36, r) = 1, N EM /I C = 10~ 3 , R = 38, Ax = 0.15, At = 5 X 10~ 4 , Q(t = 0) = 0.8, At sm = 0.1. 



ond, since the Magnus force on pinned vortices is propor- 
tional to \v — v B \ ~ r({L z )/I B — fi) at the pinning site, the 
larger radii of many of the pinned vortices in the R — 38 
system make for larger Magnus forces for smaller values of 
«£„}//. -fi). 



8 DISCUSSION 



The numerical experiments described in this paper 
demonstrate from first principles that a superfluid coupled 
to its container via vortex pinning spins down spasmodically. 
The experiments are designed to explore how the distribu- 
tions of pulsar glitch sizes and waiting times depend on pin- 
ning and stellar parameters. Computational limitations pre- 
vent us from probing the parameter regime most appropri- 
ate to the pulsar glitch problem. However, we do our best to 
mimic astrophysical conditions by adopting the same order- 
ing of dimensionless parameters as in nature. Qualitatively, 
and with respect to the canonical parameters listed in Ta- 
ble [1] the pulsar problem involves weak (but still effective) 
pinning, slow electromagnetic spin down, a huge number of 
vortices, an even larger number of pinning sites, and a heavy 
crust. 

In order to speed up the process of compiling glitch 
statistics and minimise its subjective element, we devise an 
automated glitch finder. Since the glitch-finding algorithm 
identifies all positive jumps in Afl, it is necessary to smooth 
f2 so as to filter out post-glitch oscillations or numerical 
noise. In the future, improvements in timing and observa- 
tional duty cycle should allow real pulsar spin-down curves 
to be analysed in the same way. 

Our simulations lend new insights into aspects of the 
physics of pulsar glitches. 



(1) Stronger pinning causes larger glitches, with 
(Af2g/fi) « 7.6Vi. This is not as obvious as it may sound. 
Prima facie it is equally likely that strong pinning leads to 
a large average shear with small, frequent fluctuations. In 
fact, this does not happen. 

(2) When pinning sites and vortices are equally numer- 
ous, vortices conform to the geometry of the pinning grid, 
producing larger glitches than if pinning sites are more 
abundant than vortices. Waiting times are also longer for 
fipln ~ riF than for n p i n > uf. 

(3) Glitch sizes are insensitive to increases in the pin- 



ning site abundance beyond parity (i.e. tif 



3 ), because 



vortices move across many pinning sites to negate the lo- 
cal shear, rather than stopping at an adjacent pinning site. 
The distance travelled by an unpinned vortex is a key in- 
put to heuristic glitch model s l|Warszawski fc Melatosll2008l ; 
iMelatos fc Warszaw ski 2009, for example). 

(4) A broader range of Vi results in a broader range of 
glitch sizes and more frequent glitches. 

(5) A heavier crust produces more large glitches and 
longer waiting times than does a lighter crust. 

(6) A strong spin-down torque catalyses smaller, more 
frequent glitches than a weak torque. 

(7) Pulsar astronomers derive the strength of the pulsar 
magnetic field from the spin-down rate of the crust without 
allowing for superfluid feedback. The results in this paper 
demonstrate that this is inaccurate. Even when no vortices 
have yet unpinned (e.g. t ^ 430 in the dot-dashed curve 
in Fig. [7]), Q, is a function of both the external torque and 
the superfluid dynamics [via Eq. Q]. As vortices migrate 
towards the edge of their pinning sites, the smooth, inter- 
nal spin-up torque effectively reduces the crust's spin-down 
rate. But this effect is much smaller in a real pulsar than 
in our simulations. For example, if all vortices move radially 
outward by the width of a nucleus, A(L Z )/(L Z ) w 10 -18 . 



Gross-Pitaevskii model of pulsar glitches 21 




Figure 18. Angular velocity of the crust as a function of time Q(t) for large-scale spin-down experiments with AVi/Vo = 0.0 and 0.8 
(solid and dotted curves respectively). Greyscale plots of the initial and final superfluid density |i/i| 2 for AVj/Vb = 0.0 are shown in Fig. ll7l 
Left inset: Close-up of glitch at t = 804.3 from AVi/Vo = 0.8 curve. Right inset: Close-up of glitch at t = 806.1 from AVi/V*o = 0.0 curve. 
Glitch positions are marked as arrows. Simulation parameters: Vb = 36, rj = 1, ATem/^c = 10 -3 , R = 38, Ax = 0.15, At = 5 X 10 -4 , 
n(t = 0) = 0.8, At sm = 0.1. 



quantity At am Vb AVi/Vo n pin /n F rj N EM 



(Afi/n) + + . . 

(At) + . + - 

A . + . + - + 



Table 9. Summary of how simulated glitch statistics depend on the physical levers of the model. A plus (minus) sign indicates that 
the quantity in the left column increases (decreases) when the parameter in the column heading increases. A dot indicates no clear 
correlation. 



A larger, smooth spin-up torque results from homogeneous 
outward vortex creep. 

Table [5] summarises how measures of glitch activity, 
like the mean size and waiting time and the glitch activ- 
ity parameter, depend on the underlying physical levers, 
like the pinning strength/abundance and the electromag- 
netic torque. A dot indicates where there is no/weak de- 
pendence; a plus (minus) sign indicates that the measure 
of activity increases (decreases) as the physical variable in- 
creases. Intuitively, we expect that a positive correlation 
with mean glitch size is accompanied by a positive corre- 



lation with mean waiting time. For Vb, for example, this is 
not the case; mean glitch size increases with increasing Vb, 
whereas mean waiting time does not change monotonically 
with Vb . One possible explanation is that the simulation has 
not run long enough to establish a critical velocity shear 
that ensures a net outward flux of vortices t hat matches 
the superfluid spin down to that of the crust |Alpar et alj 
11984 iMelatos fc Warszawskil [20091. The increasing vertical 
distance between the spin-down curves in Figure [7] confirms 
that this critical shear has not yet been established. 

If pinning is not strong a t the densities where only a 
neutron superfluid is present [Donati fc Pizzocherol (|2003l ) 
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Figure 19. Left: Probability density function (pdf) of fractional glitch sizes (Af2/f2) for the O(t) curves graphed in Fig. El for AVi/Vo = 
0.0 and 0.8 (solid and dotted curves respectively) unsmoothed (top) and smoothed with top-hat window function of width At Bul = 0.1 
(bottom), logarithmically binned and graphed on log-log axes. The solid grey curves are power-law least squares fits to the solid histograms 
for Af2/fi ^ 7 X 10~ 7 , with indices —1.26 and —1.23 respectively. Right: Cumulative pdf of glitch waiting times p(At) graphed on log- 
linear axes with AVi/Vb = 0.0 and 0.8 (solid and dotted curves respectively) unsmoothed (top) and smoothed with top-hat window 
function of width At sm = 0.1 (bottom). Simulation parameters: Vb = 36, r) = 1, N^m/Ic = 10~ 3 , R = 38, Ax = 0.15, At = 5 X 10 — 4 , 
fl(t = 0) = 0.8, At sm = 0.1. 



show that Vb depends on density], then glitches may arise 
from the physics of the outer pulsar core, where the neu- 
tron superfiuid coexists with a superconducting proton fluid. 
Magnetic flux tubes in the superconducting fluid compli- 
cate the superfiuid vortex dynamics, rendering the simple 
Magnus-force dynamics in this paper invalid. Hence, simu- 
lations presented in this paper are relevant only to the pulsar 
inner crust, whose density lies in the range 10 11 g cm -3 < 
P i: 10 14 S cm ~ 3 , where the major ity of protons are confined 
within nuclei (jPethick et al.ll2000l ). 

The discrete f2 jumps in the spin-down curves graphed 
throughout this paper are an alogous to the step changes in 
magnetic field seen in Fig. 1 of lField et all (|l995l ). caused by 
magnetic flux tube avalanches in type II superconductors. 
Superconductor experiments exhibit power-law size distri- 
butions over two decades, with exponents ranging from —1.4 
to —2.2. Several authors have successfully modelled the su- 
perconducting system as a self-organised critical system fluc- 
tuating about a Bean state, governed by collective unpin- 



ning and motion of flux tubes (Bass ler fc Paczuskilll998h . 
A natural theoretical framework for a self-organised critical 
system with knock-on mechanisms is a branching process. 
The characteristic po wer-law index of such systems is —1.5 
jCarreras et alj |2002l. consistently steeper than the size dis- 
tributions reported in this paper. We therefore conclude that 
either: (i) we are not simulating a self-organised critical sys- 
tem, or (ii) our definition of event size and/or identification 
algorithm is not correct. 

Although our glitch size distributions are not incon- 
sistent with power laws, we do not find conclusive evi- 
dence of collective vortex behaviour, e.g. many-generation 
cascades of un pinning events in mov i es of \ip\ 2 . This is 
not surprising. Wa rszawski fc Melatosl (|20T0bl ) found that 
sound waves from moving vortices unpin nearby vortices 
only when dissipation is low. In this paper's experiments, 
for computational exp e diency, we use 7 = 0.05, whereas 
IWarszawski fc Melatosl (|2010rJ) required 7 = 0.0025 to fa- 
cilitate an acoustic knock-on event given the same char- 
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acteristic pinning strengths as studied in this paper. In 
contrast, the proximity knock- on mechanism identified by 
IWarszawski fc Melatos! (|2010bl ) does occur in the simula- 
tions presented here. Hence collective dynamics are not ex- 
cluded completely. 

Collective effects emerge most fully in large systems. 
We therefore present preliminary results from large, compu- 
tationally expensive simulations with ~ 200 vortices. The 
large simulations yield a broader range of glitch sizes than 
smaller systems. Power-law fits to the entire size distribution 
can be rejected with high confidence. Imposing a minimum 
cut-off, to create a subset of larger glitches, yields a distribu- 
tion that can be fitted by a power law that is rejected with 
less confidence. Physically, this may be equivalent to dis- 
tinguishing between creep-like vortex motion and avalanche 
events. 

In addition to requiring larger system sizes, reliable pul- 
sar glitch simulations must also address the role of pinning 
along the length of a vortex. When a three-dimensional vor- 
tex moves radially outward, it either 'unzips' from pinning 
sites along its length or moves in small vertical segments that 
behave effectively as individual vortices. Three dimensions 
also introduce the poss ibility of turbulent superflow, i.e. a 
polarised vortex tangle |Schwarzf l988: Bar enghi et al . 2001; 
iMongiovi fc Joul 120071 ; Kobavashi fc Tsubotal |2007t ). Tan- 
gled vortices can form bundles, which behave as multiply- 
quantised vortices. A full understanding requires further in- 
vestigation using 3D simulations. 

Finally, as well as expanding GPE simulations to larger 
systems and higher resolution, we recommend that progress 
in understand ing pulsar glitches can b e achieved via cellu- 
lar automata |M orlcv & Schmidt 1996; iBassler fc Paczuskil 
1 19981 ; IWarszawski fc Melatos! 12008] ). The generic behaviour 
observed in GPE simulations can be used to devise micro- 
scopic rules for such automata. 
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